Maximally-localized Wannier functions for entangled energy bands 
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We present a method for obtaining well-localized Wannier-like functions (WFs) for energy bands 
that are attached to or mixed with other bands. The present scheme removes the limitation of the 
usual maximally-localized WFs method (N. Marzari and D. Vanderbilt, Phys. Rev. B 56, 12847 
(1997)) that the bands of interest should form an isolated group, separated by gaps from higher and 
lower bands everywhere in the Brillouin zone. An energy window encompassing A'^ bands of interest 
is specified by the user, and the algorithm then proceeds to disentangle these from the remaining 
bands inside the window by filtering out an optimally connected A-dimensional subspace. This 
is achieved by minimizing a functional that measures the subspace dispersion across the Brillouin 
zone. The maximally-localized WFs for the optimal subspace are then obtained via the algorithm 
of Marzari and Vanderbilt. The method, which functions as a postprocessing step using the output 
of conventional electronic-structure codes, is applied to the s and d bands of copper, and to the 
valence and low-lying conduction bands of silicon. For the low-lying nearly-free-electron bands of 
copper we find WFs which are centered at the tetrahedral interstitial sites, suggesting an alternative 
tight-binding parametrization. 

PACS: 71.15.-m, 71.15.Ap, 71.20.-b 



I. INTRODUCTION 

When studying electrons in solids, it is often the case 
that only a small subset of the available one-electron 
states contributes significantly to the properties under 
consideration. Moreover, the states of interest typically 
lie within a limited energy range. For instance, for mod- 
eling electron transport or magnetic properties, only the 
partially filled bands close to the Fermi energy Ep are 
needed. This is the rationale behind the tight-binding 
and Hubbard models, in which only a few energy bands 
are kept.Ela Those models rely on the existence of a min- 
imal set of spatially localized orbitals spanning the man- 
ifold of relevant states. 

In recent years there has been growing interest in ex- 
plicitly constructing such orbitals from first-principles 
density-functional calculations. One potential appli- 
cation consists in obtaining the parameters in corre- 
lated Hamiltonians by constraining the occupation of 
the orbitals to find the energy cost of deviating from 
the meaii^field solution ("constrained density- functional 
thcory"0Ll). Another arises in the context of the "dy- 
namical mean-field theory" which, when combined with 
density-functional methods, requires the specification 
of localized orbitals describing the narrow bands of 
interest .0 j-. 

Wannier functionsQ (WFs) are a very natural type 
of localized orbital for extended systems. They play a 
central role Ib, formal discussions of the tight-bindingtl 
and HubbardEl models. Traditionally they have often 
been invoked - although rarely calculated explicitly - 
as a convenient basis for describing local phenomena- 
such as impuritiesJH excitons,Ll and magnetic properties.EI 
More recently, WFs have found important applications in 



connection with lineat=scaling algorithms for electronic 
structure calculations .□ Moreover, they play an impor- 
tant role in the theory of electronic polarization and 
localization in insulators, with the former qua jitity be- 
ing related to the centers of charge 'AAg WFsHO and 
the latter to their quadratic spreads .I13E3 These develop- 
ments have also led to generalizations of the concentjOf 
Wannier functions to correlated electron systems .E^TtJ 

The main obstacles to the construction of WFs in 
practical calculations have been their nonuniqueness (or 
"gauge dependence" ) and the difficulties in dealing with 
degeneracies among the Bloch states. These have been 
overcome by the development by Marzari and Vander- 
bilt of a general and practical method for extracting 
"maximally-localized" WFs from an isolated group of 
bandsJl3 (By "isolated" we mean a group of bands that 
may become entangled with one another across the Bril- 
louin zone, but is separated from all other bands by fi- 
nite gaps throughout the entire Brillouin zone. The set 
of valence bands of an insulator constitutes an impor- 
tant example.) The method has been successfully used 
to describe the dielectric pmperties of severaLinsulating 
systems, such as crystallinell3 ancLamorphoutO semicon- 
ductors, ferroelectric pppvskites,E3 liquid wateifjij com- 
pressed solid hydrogen ,E3 and manganese oxide .Eil It has 
been implepented for plane- wa*fe|l3 linear augmented 
plane- wave and tight-bindingH3 basis sets. 

However, in many cases the group of bands of interest 
is not isolated in the above sense, especially when deal- 
ing with metals or with the empty bands of insulators. 
For example, the conduction s band of an alkali metal is 
attached at points or lines of high symmetry to higher 
bands; the d bands of a noble or transition metal are 
hybridized with an s band, which in turn is attached to 
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higher bands; the conduction bands of a copper-oxide su- 
perconductor emerge from a dense group of bands below; 
and the four low-lying antibonding bands of a tetrahedral 
semiconductor are connected to higher conduction bands. 

A successful technique that has been applied for con- 
structing localized orbitals t hat describe such bands is 
the "downfolding" techniqueE3c3 that has been devel- 
oped for electronic structure methods based on muffin-tin 
orbitals. There have also been previous attempts at con- 
structing WFs for non-isolatedpgraups of bands, namely 
for noble and tijaj^ition metalaHjE] and for tetrahedral 
semiconductors B3'E2l These attempts fall into two cate- 
gories: (i) the WFs are obtained directl*i-|from a vari- 
ational principle, as suggested by Kohn,E3 or (ii) they 
are obtained as Fourier transforms of Bloch functions, 
with the help of a model Hamiltonian that reproduces 
the band structuce in the desired energy range, as sug- 
gested by Bross.Eil 

We will describe an alternative Wannier-based ap- 
proach that is closer in spirit to the Fourier transform 
method of Bross and co-workers, but does not require 
the construction of an auxiliary model Hamiltonian. The 
method can be regarded as an extension to the case of 
attached bands of the mappially-localized WF method 
of Marzari and Vanderbilt.ll3 It has the desirable features 
that it can be implemented with any basis set (e.g., plane 
waves), and requires minimal user-intervention (the only 
"adjustable parameter" being a specification of the en- 
ergy range of interest). Like the approach of Ref. |l^, 
ours is a "postprocessing" method, taking as its input the 
Bloch eigenstates and eigenvalues calculated by a stan- 
dard electronic-structure code. 

Strictly speaking, the resulting orbitals are not WFs 
(or even "generalized WFs"ll3) in the usual sense. They 
are nevertheless Wannier-like in the fundamental sense 
that they are obtained via an integral over the Brillouin 
zone of Bloch-like functions. As such they form an or- 
thonormal, localized basis of the same Bloch subspace 
from which they were constructed. 

The power of the present approach is illustrated by one 
particularly striking result that emerged from the work. 
In Sec. IV B 3 we find that a rather natural representa- 



tion of the low- lying bands of an fee metal like copper can 
be made in terms of a set of five Cu c?-like WFs and two 
additional WFs centered at the tetrahedral interstitial lo- 
cations. This provides a basis for a concise tight-binding 
representation of copper that has not, to our knowledge, 
previously been considered. 

The paper is organized as follows. In Sec. |l| we review 
the method of Marzari and Vanderbilt for obtaining well- 
localized WFs for an isolated group of bands. In Sec. Ill 
we describe our procedure for dealing with attached en- 
ergy bands, and in Sec. tV we illustrate it with a set of 
applications. Finally, in Sec. ^ we present a summary 
and conclusions. 



II. MAXIMALLY-LOCALIZED WANNIER 
FUNCTIONS FOR AN ISOLATED GROUP OF 
BANDS 

A set of WFs ?i'„R(r) — Wn{r — R) labeled by Bra- 
vais lattice vectors R can be constructed from the Bloch 
eigenstates ^„k of band n using the unitary transforma- 
tion 
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where v is the volume of the unit cell of the crystal and 
the integral is over the Brillouin zone. Except for the con- 
straint ^pn,k+G = '<Pnk for all reciprocal lattice vectors G, 
the overall phases of the Bloch functions V'nk = e*'^ '"u„k 
are at our disposal. However, a different choice of phases 
(or "gauge"). 
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does not translate into a simple change of the overall 
phases of the WFs; their shape and spatial extent will in 
general be affected. If the band is isolated, Eq. (||) is the 
only allowed type of gauge transformation for changing 
the set of WFs w„ (r — R) associated with that band. In 
the case of an isolated group of N bands, the allowed 
transformations are of the more general form 
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where U^^"^ is a unitary matrix that mixes the bands at 
wave vector k. The resulting orbitals are called "gener- 
alized Wannier functions" Jl3 

Once a measure of localization has been chosen and an 
isolated group of bands specified, the search for the cor- 
responding set of "maximally-localized" WFs becomes a 
problem of functional minimization in the space of the 
matrices U^^l The strategy of Ref. |l| consists in mini- 
mizing the sum of the quadratic spreads of the Wannier 
probability distributions |?«„(r)|^. 
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where the sum is over the chosen group of bands and 
(^)n ~ I ^\u}n(jc)\'^ dr, etc. Interestingly, the resulting 
"maximally- localized" (or "maxloc") WFs turn out to 
be real, apart from an arbitrary overall phase factor. 

In numerical calculations the Bloch states V'nk are com- 
puted on a regular mesh of fc-points in the Brillouin zone; 
the integral in Eq. ([T]) is then replaced by a sum over the 
points in the mesh. In Ref. |l^ an expression was derived 
for the gradient of the spread functional with respect 
to an infinitesimal rotation SU'^^^ of the set of Bloch or- 
bitals. The only information needed for calculating the 
gradient are the overlaps 
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where b are vectors connecting a mesh point to its near 
neighbors. Once the gradient is computed, the mini- 
mization can proceed via a steepest-descent or conjugate- 
gradients algorithm. 

In Ref. the spread Q was decomposed into two 
terms, 



(6) 



both of them non-negative. The first measures the fc- 
space dispersion of the band projection operator, while 
the second reflects the extent to which the Wannier func- 
tions fail to be eigenfunctions of the band-projected po- 
sition operators, fli will play a central role in the present 
work. For an isolated group of bands it is invariant un- 
der any gauge transformation so that minimizing fl 
amounts to minimizing fl. When using a regular mesh of 
fc-points, r^i is given by 
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where N\^p is the total number of /c-points, N is the num- 
ber of bands in the group, and wi, is a weight that arises 
from the discretization procedure by which derivatives, 
with respect to k are approximated by finite differences. t£l 
The corresponding expression for Q can be found in 
Ref. 



III. MAXIMALLY-LOCALIZED WANNIER 
FUNCTIONS FOR ATTACHED BANDS 

A. Description of the method 

For definiteness let us suppose we want to "disentan- 
gle" the five d bands of copper from the s band which 
crosses them (see Fig. 0) and construct a set of well- 
localized WFs associated with the resulting d bands. 
Heuristically the d bands are the five narrow bands and 
the s band is the wide band. The difficulty arises because 
there are regions of A:-space where all six bands are close 
together, so that as a result of hybridization "the distinc- 
tion between d-band and s-band levels is not meaningful" 
(Ref. I, p. 288). 

Let us now outline our strategy, which can be divided 
in two steps. First we cut out an energy window that en- 
compasses the TV bands of interest {N = 5 in our exam- 
ple). Figs. |l|(a) and |l|(b) correspond to different choices 
for this energy window. At each /c-point the number 
of bands that fall inside the window is equal to or larger 
than the target number of bands N. This procedure de- 
fines an A^k-dimensional Hilbert space .F(k) spanned by 
the states it„k within the window. If at some k = N, 
there is nothing to do there; if A^k > N our aim is to find 



> 



c 



■ (a) 


y 






\, 























> 



c 



- (b) 




\ 












\ 
\ 


/ 



r X w L r K 

FIG. 1. Solid line: Calculated band structure of cop- 
per. Dotted line: Interpolated bands obtained from the five 
d-like Wannier functions, (a) and (b) differ in the choice of 
the energy window used to compute the Wannier functions 
([-9.59,-0.29] eV in (a) and [-9.59,7.21] eV in (b)). The 
zero of the energy scale is at the Fermi energy. 

the A^-dimensional subspace iS(k) C J-'(k) that, among 
all possible A^-dimensional subspaces of T{^^ leads to 
the smallest 57i (Eq. (0)). (Recall that for an isolated 
group of bands rJi is gauge-invariant, since it is an intrin- 
sic property of the manifold of states. Thus fii can be 
regarded as a functional of 5(k).) In the second step we 
work within the optimal A^-dimensional subspaces iS(k) 
selected in the first step, and minimize 17 using the al- 
gorithm of Marzari and Vanderbilttj summarized in the 
previous section. The end result is a set of N maximally 
localized WFs and the corresponding N energy bands. 
We emphasize that it is the first step (minimization of 
i7i) that is new with respect to Ref. Qq. 



B. Physical interpretation of fJi 

Why is minimizing f2i a sensible strategy for picking 
out the d-bands? This can be understood by noting that 
heuristically fii measures the "chaHge of character" of the 
states across the Brillouin zoneJlj Indeed, Eqs. (|^) and 
(^ show that ill is small whenever | (unk|um.k+b) the 
square of the magnitude of the overlap between states 
at nearby fe-points, is large. Thus by minimizing fii 
we are choosing self-consistently at every k the subspace 
5(k) that has minimum "spillage" or mismatch (see be- 
low) as k is varied. In the present example this optimal 
"global smoothness of connection" will be achieved by 
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keeping the five well-localized d-like states and exclud- 
ing the more delocalized s-like state. We will gain more 
intuition about the meaning of minimizing fli while dis- 
cussing specific examples in Sec_-|M 

What is meant by "spillage" t°fE3 becomes clear once 
we rewrite Eq. (Q) as 



k„ A 



with 
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where Pk = J2n l'"nk)(Mnk| is the projector onto 5(k), 
Qk = and the band indices to, n run over 1, . . . , N. 

7k, b is called the "spillage" between the spaces 5(k) and 
5(k -I- b) because it measures the degree of mismatch 
between them, vanishing when they are identical. 

Further discussion of the geometrical and physical in- 
terpretation of ill can be found in Refs. ^ and ^ In 
particular, it has been shown that the value of Qi asso- 
ciated with the valence bands of an insulator is the ex- 
perimentally measurable mean-square quantum fluctua- 
tion of the ground state macroscopic polarization. □ This 
can be interpreted as the quadratic spread of an appro- 
priately defined collective center-of-mass distribution for 
the electrons, and can be recast as an electronic localiza- 
tion length squared. Hence our procedure of minimizing 
ill selects the A^-dimensional subspaces 5(k) where the 
electrons are most localized in the above sense (assum- 
ing for the purpose of this argument that all the electron 
states in those subspaces are occupied). 

Finally we note in passing that our two-step procedure 
of minimizing first f2i and then n is in principle different 
from directly minimizing their sum fl. In view of the dis- 
cussion presented above, we believe that the procedure 
adopted here is conceptually the more natural of the two, 
although we would expect them to yield similar results 
in practice. Also, as we will now show, the separate min- 
imization of f2i turns out to be a particularly simple and 
robust procedure. 



C. Iterative minimization of Qj 

Since the functional that we wish to minimize cou- 
ples states at different fc-points, the problem has to be 
solved self-consistently throughout the Brillouin zone. 
Our strategy is to proceed iteratively until the optimal 
"global smoothness of connection" is achieved. On the i- 
th iteration we go through all the fc-points in the grid, and 

(i) 

for each of them we find N orthonormal states u^j^, defin- 
ing a subspace iS*^*''(k) C jF(k) such that the "spillage" 
over the neighboring subspaces iS'*~^''(k + b) from the 
previous iteration is as small as possible (Fig. 0). 
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FIG. 2. Schematic representation of the subspaces of 
Bloch-like states on a grid of fc-points. Our procedure consists 
of iteratively minimizing the "spillage" , or degree of mismatch 
(see text), between the subspaces at neighboring fc-points. 



Using Lagrange multipliers to enforce orthonormality, 
the stationarity condition at the i-th iteration reads 
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where A^*' is an x iV matrix. Let 
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where, according to Eq. (||), 

a.«(k) = 5:u;,r« 
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The first term in Eq. (10) now becomes 
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From Eq. (|T^) we find 
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where f'k^ J'' is the projector onto 5'' ^-'(k -I- b). Like- 
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wise, one easily obtains 
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Combi ning the previous equations, the stationarity con- 
dition ( |lO[) becomes 



N 



where A-l^^.]^ = i^^p/'^)^nL,k- choosing a unitary 
transformation that diagonahzcs A^^'' , this can be recast 
as an eigenvalue equation: 
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The eigenvalues of the above equation obey < X^^y. < 
J2h ^b', in particular, A^^j^ < Ylib ^fc whenever the eigen- 
state mI^L does not lie completely within all of the nearby 



subspaces 5^* ^^(k + b). Combining Eqs. ( p^ ) and (17) 
we find 



AT 



C.«(k) 



^ E - E ^ 



mk" 



(18) 



It is clear from Eqs. ( [Ill ) and ( |l8| ) that when construct- 
ing 5^'-* (k) one should pick the N eigenvectors of Eq. ( p7| ) 
with largest eigenvalues, so as to ensure that the station- 
ary point corresponds to the absolute minimum of 

Self-consistency is achieved when 5^'-'(k) = iS'-*~^-'(k) 
at all the grid points. We have encountered cases where 
the iterative procedure outlined above was not stable. 
In those cases, the problem was solved by using as the 
input for the present step a linear mixing of the input and 
output subspaces from the previous step. More precisely, 
the eigenvalue equation (jlj) was replaced by 
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with A typical value is Q!=0.5. 

In practice we solve Eq. (|l^) in the basis of the origi- 
nal A^k Bloch eigenstates M„k inside the energy window. 
Each iteration then amounts to diagonalizing the follow- 
ing A'k X A^k Hermitian matrix at every k: 
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Since these are small matrices, each step of the itera- 
tive procedure is computationally cheap. In particular, 
the time-consuming computation of the overlap matri- 
ces M^'^'^^ of Eq. (H) can be done once and for all at 
the beginning, using the original Bloch eigenstates in- 
side the energy window; their subsequent update during 
the iterative minimization is very inexpensive. An anal- 
ogous situation occurs when updating the matrices U^^^ 
in Eq. (^ duripg the minimization of f2 to obtain the 
"maxloc" WFs.EJ 



D. Initial guess for the subspaces 

In order to start the iterative minimization of f2i, the 
user should provide an initial guess for the subspaces 
S (k) . We have found that the minimization procedure is 
quite robust, in the sense that it is able to arrive at the 
global minimum starting from a very rough initial guess. 
In practice we usually select the initial subspaces follow- 
ing a strategy very similar to the one outlined in Ref. ^ 
for starting the minimization of ^l. 

A set of A^ localized trial orbitals gra(r) is chosen cor- 
responding to some rough initial guess at the WFs, and 
these are then projected onto the A'k Bloch eigenstates 
inside the energy window. 



yjik, 



= E^^ 



mk/ 



(22) 



where Amn = {"iprnk |<?n) is an A^k x N matrix. The re- 
sulting A^ orbitals are then orthonormalized via Lowdin's 
symmetric orthogonalization procedure,tj i.e.. 
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where S„in = (0mk|(Ank) = (^^^)™„- Finally these 
Bloch-like functions are converted to cell-periodic func- 
tions = ^~*'""^ik- matrix AS~^^'^ can easily 
be computed by pecforming the singular-value decompo- 
sition A — ZDV^ where Z and V are A'^k x A'^k and 
N X N unitary matrices respectively, and D is A'k x N 
and diagonal. This leads to AS~^^'^ = ZIV, where 1 is 
the Ak X A^ identity matrix. 



E. Minimization of SI 

At the end of the first step of our procedure (mini- 
mization of f2i) we are left at each A:-point with an A^- 
dimensional subspace 5(k), and for definiteness we diag- 
onalize the Hamiltonian inside this subspace to obtain A^ 



Bloch-like eigenfunctions tjjnk 



o^k-r 



■''u„k and eigenvalues 



e„k- The second step is to find the N x N unitary ma- 
trices U^^^ (Eq. (||)) that, applied to the V'nk, produce 
the rotated set of Bloch-like states that is transformed 
via (|l|) into the maximally-localized WFs w„r. X14s is 
done using the method of Marzari and VanderbiltEEl for 
minimizing f2, briefly discussed in Sec. ||. An initial guess 
for the unitary matrices J/'-'^-' is obtained by projecting a 
set of A^ localized orbitals onto the states tpnk- Typically 
the same set of orbitals is used as in the initialization step 
for the minimization of Qi. (In our experience, when a 
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particularly bad choice of trial orbitals is made, the min- 
imization of fli is less likely to become trapped in local 
minima than the minimization of f2.) 

F. Interpolated band structure 

Starting from the "maxloc" WFs, the corresponding 
energy bands can be computed at arbitrary points in 
the Brill ojii.p .zone using a Slater-Koster interpolation 
scheme. EIlLiJcj Of course, the interpolation could proceed 
directly from the non-rotated states u„k; however, use of 
the optimally rotated ones ensures that -tlie interpolated 
band structure is as smooth as possible.L3 

The interpolation procedure involves first calculating 
the Hamiltonian matrix for the rotated states, 



(24) 



where i7m„(k) = 

Emk^m.n- Ncxt wc Fouricr transform 
^(rot)(li) into a set of A^kp Bravais lattice vectors R 
within a Wigner-Seitz supercell centered around R = 0: 



(25) 



k 

= {Wmo\H\wnR), 



where H is the effective one-particle Hamiltonian. Fi- 
nally we Fourier transform back to an arbitrary fc-point, 



R 



(26) 



and diagonalize the resulting matrix to find the interpo- 
lated energy eigenvalues. 



G. Inner energy window 

In some situations one wants to construct orbitals that 
describe the original bands exactly only in a limited en- 
ergy range. This can occur when studying transport 
properties for which only the states within some small 
energy range of the Fermi level (say, ±1 eV) are relevant. 
The challenge is to construct orbitals that achieve that 
goal while remaining as localized as possible. What the 
resulting interpolated bands look like outside the energy 
range of interest is largely immaterial, since it will not af- 
fect the low-energy physics. (Typically they will tencLtp 
remain close in energy to the target range of interest .l3) 

A simple extension of the formalism described in the 
previous sections can produce such orbitals. The idea is 
to introduce a second ("inner") energy window - con- 
tained within our original ("outer") window - inside 
which the original bands are to be described exactly. Let 
Mk be the number of bands that fall within the inner win- 
dow at k, so that Mk < N < TVk. Then we have to mini- 
mize fii under the constraint that the Mk original Bloch 



states inside the inner window must be included in the 
subspace iS(k). We are therefore only free to choose the 
remaining N — Mk states when constructing iS(k) . Those 
will have to be extracted from the subspace spanned by 
the A'k — Mk original Bloch cigenstates that are inside 
the outer window but outside the inner window. That 
can be achieved by a straightforward modif ication of the 
iterative procedure described in Sec. [ill C 

(k) in Eq. ^ becomes an (iVk 
matrix, and we pick the N — A/k leading eigenvectors 

The only remaining iss ue is how to modify the initial 
ization procedure of Sec 



The matrix 
A/k) X (iVk - Mk) 
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in order to accommodate 
the inner window. Since the first Mk basis vectors of 
the trial subspaces S (k) are predetermined, we want the 
modified procedure to provide the remaining N — M]^ vec- 
tors. Let Q(k) be an A^-dimensional space obtained by 
projecting the N trial orbitals onto the A^k s tates inside 
the outer window, as described in Sec. [Ill D| . Let Pg;(k) 
be the A^k x -^k matrix that is the projection operator 
onto ^(k) as expressed in the space .7-"(k). Similarly, 
define Pi„nor(k) as the A'k x A'k projection matrix onto 
the inner window states, and (5inncr(k) = 1 — Pinner (k). 
Then choose the remaining N — Mk basis vectors to be 
the eigenvectors corresponding to the N — Mk largest 
eigenvalues of 



(k)Pg(k)(5inncr(k)|w) — X\v) . 



(27) 



Such vectors have the desired properties: (i) They are 
orthogonal to the states inside the inner window, and 
(ii) because A — (w|Pg;(k)|w), it is clear that by choosing 
the eigenvectors with the largest eigenvalues we guaran- 
tee that their overlap with the space ^y(k) is as large as 
possible, while satisfying the constraint (i). 

Other kinds of constraints on the minimization of fij 
may also be useful. For instance, one might want to "pin 
down" the desired bands at high-symmetry fc-points to 
ensure that the interpolated bands coincide with them at 
those points. 



IV. RESULTS 

A. Computational details 

The calculations were performed within the local- 
density approximation to density-functional theory, us- 
ing a plane-wave basis set and TrouUier-Martins norm- 
conserving pseudopotentiala^ in the Kleinman-Bylander 
representation. The energy cutoff was set to 75 Ry for 
copper and 35 Ry for silicon, and the lattice constants 
were 6.822 bohr and 10.260 bohr respectively. The com- 
puted self-consistent Bloch eigenfunctions and eigenval- 
ues that fell inside the prescribed energy window were 
stored to disk. They were used as the input for the min- 
imization of r^i, which was carried out as a separate. 
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(a) 



(b) 




FIG. 3. Contour-surface plots of the two Eg Wannier func- 
tions associated with the "disentangled" d bands of copper 
shown in Fig. ^(b). The amplitudes are +0.5/ ^/v (light gray) 
and —0.5/^11 (dark gray), where v is the volume of the prim- 
itive cell. 

postprocessing operation. This produced an optimal sub- 
space characterized by a new set of N Bloch eigenfunc- 
tions and eigenvalues per fc-point, which were taken as 
the input for constructing the "maxloc" WFs and the 
interpolated bands. In all the cases we have found the 
"maxloc" WFs to be real (apart from an overall phase 
factor), as was alreadi£-|the case when dealing with iso- 
lated groups of bandsJla The self-consistent calculations 
were performed on a 10 x 10 x 10 Monkhorst-Pack mesh 
of fc-points for copper, and 6x6x6 for silicon. During 
the minimization of fli and a 10 x 10 x 10 uniform 
grid was used for both copper and silicon. This grid was 
shifted in order to include the F point (fc — 0), so as 
to ensure that the "maxloc" WFs have the desired sym- 
metry properties among themselves. (For instance, if a 
grid is used for silicon that does not include F, the four 
antibonding WFs in a unit cell do not all have the same 
spread.) The mixing parameter a in Eq. ( pO| ) was set to 
0.5. 



B. Copper 

Wannier functions for noble and transition metalSriaKe 
previously been computed using various approaches. c^EII 



Below, taking copper as an example, we show how the 
present scheme can be used to "disentangle" the narrow 
d bands from the nearly- free-electron bands, allowing us 
to treat each group of WFs separately. Alternatively, 
one can also treat the narrow and the nearly-free-electron 
bands as a single group. 

1. Narrow d bands 

First, an energy window was chosen such that at each 
fc-point in the grid it contained six or seven energy eigen- 
values. As indicated in Fig. |l|, the precise range of the 
window is largely at our disposal; unless explicitly stated 
otherwise, the numbers given below pertain to Fig. |l|(b). 
In order to extract the five d bands, we set TV = 5 and 
initialized the minimization of both fli and fl from five 
trial Gaussians of r.m.s. width 1 bohr, each modulated 
by a different I = 2 angular eigenfunction. After ~ 50 
iterative steps fij was fully converged, having decreased 
from an initial value of 9.957 bohr^ to 8.483 bohr^. Dur- 
ing the subsequent minimization of Q the total Wannier 
spread decreased only slightly, from 8.563 bohr^ to 
8.556 bohr^. In agreemeajc with previous experience on 
isolated groups of bands ,t3 we found for the d bands that 
at the minimum fii 3> f^. 

The bands obtained by interpolation using the five 
"maxloc" WFs are shown as dotted lines in Fig. |l|, to- 
gether with the original band structure. As expected, 
whenever the dispersive s-like band is far from the narrow 
d bands, so that they retain their separate identities, the 
interpolated bands are very close to the narrow bands. 
However, whenever the six bands are close together, and 
thus strongly hybridized, the interpolated bands remain 
narrow, which suggests that they are mainly d-like in 
character. (Heuristically they can be viewed as the bands 
obtained by artificially "switching off" the Hamiltonian 
matrix elements between s and d WFs, i.e., by remov- 
ing the hybridization.) The d character is confirmed by 
inspection of the contour-surface plots of the "maxloc" 
WFs, two of which are shown in Fig. ^. The quadratic 
spreads of the five WFs are not exactly equal, because 
of the Cg — t2g splitting of the d-states; those shown in 
Fig. H {eg orbitals) have a spread of 1.700 bohr^ each, 
whereas the remaining three {t2g orbitals) each have a 
spread of 1.718 bohr^. These numbers are only slightly 



TABLE I. Variation of the optimal Wannier spread SI and 
its gauge- invariant part f2i (in bohr^) with the choice of en- 
ergy window range (in eV), for the d bands of copper. 



Window ranf 


;e 


Total spread 




Min 


Max 


Qi 


n 


-9.59 


-0.29 


15.373 


16.489 


-9.59 


2.21 


10.404 


10.621 


-9.59 


7.21 


8.483 


8.556 


-9.59 


12.21 


7.634 


7.667 
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FIG. 4. (a) Dotted lines: the s-d bands of copper obtained 
by extracting the optimal six-dimensional subspace 56(k) in- 
side the window, (b) Dotted lines: d bands associated with 
optimal five-dimensional subspace 55(k) C 56(k). Dashed 
line: s band 5i(k) isolated by taking the complement of 
55 (k). 



larger than the ones reported in Table III of Ref. |2^, ob- 
tained using a different method and a sparser sampling 
of the Brillouin zone. 

In our procedure there is one adjustable parameter, 
namely the range of the energy window. This range 
should be wide enough that it encompasses the bands 
of interest, but not be so wide that it also includes other 
bands of similar character (e.g., higher d bands). In the 
limit of a very wide window the spaces ^(k) would con- 
tain a complete set of states, so that by mixing in states 
far away from the energy range of interest but of similar 
character, the spread of the WFs could be made arbi- 
trarily small (and the corresponding bands would become 
flat). Table | shows how the optimal Wannier spreads are 
affected by varying the window range within reasonable 
bounds. As anticipated, the spread decreases with in- 
creasing energy rangeEj The change in the interpolated 
energy bands is less pronounced, although they do be- 
come somewhat narrower (compare Figs. |^(a) and|^(b)). 
In particular, the upward shift of the lowest interpolated 
band at L is caused by mixing with the^aeventh band, 
which has the same symmetry label (Li).c3 



2. Nearly- free- electron hand 

The unconstrained minimization of fJi usually pro- 
duces narrow bands, since the character of the Bloch 
states in such bands tends to have only a small vari- 
ation across the Brillouin zone, corresponding to well- 



TABLE II. Spreads of the "maxloc" WFs for the separate 
d-band and s-band subspaces (S^ and Si), and for the com- 
bined s-d subspace S&. The numbers in parentheses are the 
f2i values, and t stands for tetrahedral-interstitial-centered 
orbital. The corresponding bands are displayed in Fig. W. 



Two separate subspaces 



One combined subspace 





1.710 




deg 


1.731 




1.710 




deg 


1.731 


dt2g 


1.808 




dt2g 


2.328 


dt2g 


1.808 






2.328 


dt2g 


1.808 






2.254 




8.844 


(8.745) 






t 


12.929 




t 


10.263 


^^min[5l] 


12.929 


(10.826) 




20.634 (16.506) 



localized electrons (this may not be the case in the pres- 
ence of avoided crossings). The method is therefore ide- 
ally suited for directly extracting the narrow d bands 
from the s-d complex. If instead one is interested in iso- 
lating the wider, nearly-free-electron s band, direct min- 
imization of fli for one-dimensional subspaces is not the 
appropriate strategy. Instead one can proceed as follows. 
First choose an energy window that includes the s-d band 
complex (we used the one indicated in Fig. |](b)). Then 
minimize fii with A'^ = 6; this produces a six-dimensional 
subspace iSeC*) throughout the Brillouin zone that con- 
sists of the s-d band complex. Next extract the five d 
bands by minimizing Qi within iS6(k) choosing N — 5; 
this yields a space ^5 (k) C (k) . The difference between 
the two is a one-dimensional space Si (k) containing the 
desired band. Fig. ^(a) shows the bands associated with 
56(k), and Fig. H(b) shows the bands corresponding to 
S5{k) and 5i(k). 

In Table ^ are presented the optimal Wannier spreads 
for the different subspaces. We find that the spread of the 
s-like WF is considerably smaller than the ~45 bohr^ re- 
ported in Table III of Ref. |2^. Moreover, contrary to what 
one might have expected, that WF is centered not on an 
atom, but on a tetrahedral interstitial site, as shown in 
Fig. g|(a). Since there are two such sites per atom, a 
breaking of symmetry must have occurred when select- 
ing the subspace 56(k). Indeed there are two degenerate 
minima of Qi with = 6, one for each of the interstitial 
sites. If the minimization is initialized by projecting five 
d-like orbitals plus one s-like orbital, all atom-centered, 
the breaking of symmetry occurs spontaneously during 
the iterative procedure (the minimization of f2i reaches 
a plateau, presumably a saddle point, and eventually the 
algorithm finds its way towards one of the two minima). 
If instead the s trial orbital is centered around one of 
the tetrahedral interstitial sites, the minimization starts 
inside the basin of attraction of the corresponding mini- 
mum. 

Finally, as a sim ple illustration of the "inner window" 
idea of Sec. 



Ill G, we show in Fig. H the single band 



(A^=l) that results when an inner window is selected in 
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(c) 




FIG. 5. Contour-surface plots of interstitial-centered 
"maxloc" WFs. (a) t-like WF associated with the subspace 
56 (k) of Fig. I and Table [l[ (b) WF associated with the band 
in Fig. ^; (c) t-like WF associated with the subspace 57(k) in 
Fig. 0(a) and Table ^ The amplitudes are -1-0.5/^ (light 
gray) and —0.17/^, —0.3/y/v, and — 0.25/v^ (dark gray) in 
(a), (b), and (c) respectively. 



r X w L r K 

FIG. 7. (a) Dotted lines: Interpolated bands associated 
with the optimal subspace 57(k) containing five d-like WFs 
and two tetrahedral-interstitial-centered WFs. (b) Dark dot- 
ted lines: d bands associated with optimal five-dimensional 
subspace 55(k) C 57(k). Light dotted lines: dispersive bands 
52 (k) isolated by taking the complement of 55 (k). 

the energy range below the d bands. As expected, the 
interpolated band is identical to the original one inside 
that window. Moreover, it remains quite narrow outside, 
where it acquires a pronounced d character. (This means 
that the cost in fij of changing from s to ci character is 
more than compensated by the smaller dispersion - and 
hence smaller ili - of the more localized d-like states.) 
Accordingly, the "maxloc" WF, shown in Fig. ||(b), is 
again centered at a tetrahedral interstitial site, like the 
WF of Fig. |( a), but now it has a substantial admixture 
of d-like satellites and a smaller spread, fl = 7.323 bohr^ 
{ni = 7.306 bohr^). 

The results of this Section indicate that the occurrence 
of a symmetry breaking in the minimization of Qi with 
a "maxloc" WF centered at a tetrahedral interstitial site 
appears to be a rather robust result. 



> 



5 










\ 
> 





_ 








/ 

/ 


5 








\ 

— A 

















r X w L r K 

FIG. 6. Dashed line: Band obtained using both an inner 
and an outer energy window. 



3. Symmetric two- WF description of dispersive bands 

Remarkably, we find that the symmetry can be re- 
stored, and a more faithful overall description of the 
bands can be achieved, by bringing in just one more dis- 
persive band and working with a set of seven WFs. More 
precisely, we choose an energy window such as the one in- 
dicated in Fig. 0(a), containing seven or more bands, and 
minimize fii with N — 7. (To ensure that the low-energy 
part of the band complex is well described, we freeze it 
inside an inner window.) After applying the localiza- 
tion procedure we obtain, besides the five d orbitals, two 
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TABLE III. Spreads of the "maxloc" WFs for the sep- 
arate d-band and low-lying dispersive bands subspaces (55 
and 52), and for the combined subspace Sr- The numbers 
in parentheses are the Qi values, and t stands for tetrahe- 
dral-interstitial-centered orbital. The corresponding bands 
are displayed in Fig. m 



Two separate subspaces One combined subspace 





1.687 




deg 


1.687 


deg 


1.686 




deg 


1.687 


dt2g 


1.472 




dt2g 


1.737 


dt2g 


1.472 




dt2g 


1.737 


dt2g ^ 


1.472 




dt2g 


1.737 


^min ['^5] 


7.788 


(7.751) 






t 


8.568 




t 


7.812 


t 


8.568 




t 


7.812 


Slmin [52] 


17.136 


(16.822) 


nmin[57] 


24.209 (22.034) 



equivalent WFs, each centered at one of the two tetra- 
hedral interstitial sites. One of the latter is shown in 
Fig. Jq(c) . The optimal Wannier spreads are given in Ta- 
ble ffl; it can be seen that the spread of each of the two 
interstitial WFs is considerably smaller than that of the 
single interstitial WF in Table || and Fig. ||(a). 

Fig. 0(b) shows the d-like bands associated with the op- 
timal five-dimensional subspace 'S5(k) C S-j{k), as well as 
the dispersive bands associated with 52(k), the comple- 
ment of 55(k) inside iS7(k). There is an upward shift in 
energy of the states X3, W3, and Li in the narrow bands, 
due to mixing with the states of the same symmetry in 
the dispersive bands, which suffer a downward shift of 
the same magnitude. 

The fact that our procedure naturally generates a pair 
of WFs centered at the tetrahedral interstitial sites can 
be rationalized in terms of a tight-binding description of 
the nearly-free electron states. The tetrahedral intersti- 
tial sites form a simple cubic lattice, so that in view of 
Fig. ^(c) one might imagine that the electronic states 
of these WFs would be roughly analogous to those of a 
nearest- neighbor tight-binding model of s orbitals on the 
sites of a simple cubic lattice. Indeed we have checked 
that the main qualitative features of the interpolated 
bands associated with the two interstitial-centered WFs 
(light dotted lines in Fig. 0(b)) are captured by such a 
tight-binding model, but folded back into the fee Bril- 
louin zone to give two bands instead of one. 

The quality of the interpolated bands in Fig. 0(a) 
suggests that the two tetrahedral-interstitial-centered or- 
bitals (which we denote as t orbitals) complement the five 
atom-based d orbitals nicely to form a basis {t^d^) for a 
tight-binding parametrization of the copper bands. This 
requires only oneanore basis function than the traditional 
"minimal basis"cd sd^ (five d plus one s atomic orbitals), 
while|-S|till remaining more economical than the sp^d^ 
basis.E3 The three bases are compared in Table At 
each high-symmetry fc-point we list, in order of increasing 
energy, the symmetry labels of the states that occur in 



TABLE IV. A list, in order of increasing energy, of the 
symmetry labels of selected states in the band structure of 
copper (taken from Ref. and whether or not they are 
captured by each of the tight-binding bases discussed in the 
text. An asterisk (*) indicates that the state is occupied. 





Degeneracy 




t^d:' 


sp^d^ 


Fi 


1* 


yes 


yes 


yes 


F25' 


3* 


yes 


yes 


yes 


F12 


2* 


yes 


yes 


yes 


F2' 


1 


- 


yes 


- 


Fl5 


3 


— 




yes 


Xi 


1* 


yes 


■ 


yes 


X3 


1* 


yes 


yes 


yes 


X2 


1* 


yes 


yes 


yes 


Xs 


2* 


yes 


yes 


yes 


X,, 


1 


- 


yes 


yes 


Xi 


1 


yes 


— 


yes 


Xs' 


2 


— 


— 


yes 


X-A 


1 


— 


yes 




Li 


1* 


yes 


yes 


■ 


L3 


2* 


yes 


yes 


yes 


L3 


2* 


yes 


yes 


yes 


L2' 


1* 


— 


yes 


yes 


Li 


1 


yes 


yes 


yes 




1 








Ly 


2 






yes 


W2' 


1* 


yes 


yes 


yes 


W3 


2* 


yes 


yes 


yes 


Wi 


1* 


yes 


yes 


yes 


Wv 


1* 


yes 


yes 


yes 


W3 


2 




yes 


yes 


Wv 


1 






yes 


Wi 


1 


yes 




yes 



a detailed band-structure calculation (e.g., Ref. ^), and 
then whether or not they are captured by each of the 
tight-binding bases. Inspection of the table clarifies that 
the t^d^ basis has some very attractiszje features. Whereas 
the sd^ basis misses the X^i statald (unoccupied p-like 
state not far above Ep) and, even more importantly, the 
L2' state (occupied p-like state just below Ep), t^d^ gets 
the symmetries right up to at least the first state above 
Ep at each high-symmetry fc-point. Even sp^d^ does 
not do this, failing at the F point, since the state T2' 
has / character. A consequence of this analysis is that 
the t orbitals cannot be constructed solely from s and 
p orbitals. This can also be seen from Fig. ||(c): The 
positive-amplitude central portion of the WF can be in- 
terpreted in terms of a superposition of four sp hybrids 
coming from each of the four surrounding copper atoms 
and pointing towards the interstitial; however this pic- 
ture cannot account for the six negative lobes. 

To conclude, we note that the sp^d^ description can 
also be obtained from our procedure, by minimizing fli 
with N — 9 within a window containing eleven or more 
bands (e.g., with the upper bound at 32.2 eV). The 
"maxloc" WFs are then five atom-centered d-like orbitals 
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plus four equivalent sp^-like hybrids centered near the 
atom. 



C. Silicon 

Several authors have previously discussed and com- 
puted WFs for silicon and other tetrahedral semiconduc- 
tors. Some works have|_£acnsfid_pn the WFs associated 
with the valence bands Jlj'cElacfl while othep,iL*ve also 
dealt with the lowest four conduction bands.P'-' 



1. Bond orbitals 

A set of eight bond-centered WFs, four bonding and 
four antibonding, can be obtained by using separate en- 
ergy windows for each of the two groups, as indicated 
in Fig. ^(a). Since the valence bands form an isolated 
group, inside the corresponding window A^k = iV = 4 
throughout the Brillouin zone. Hence there is no free- 
dom for minimizing Oj, and one can proceed directly 
with the minimization of fl to compute the "maxloc" 
WFs, as done in Ref. The resulting bands are es- 
sentially indistinguishable from the original ones, since 
for such a dense fc-mesh the interpolation error is very 
small. The trial orbitals used to start the minimization 
were bond-centered Gaussians with a root mean-square 
(r.m.s.) width of 1.89 bohr. The value of the optimal 
spread was ft = 30.13 bohr^, of which 28.39 bohr^ came 
from r^i. 

The use of an energy window becomes necessary for 
the four low- lying empty bands, which are attached to 
higher bands. As trial orbitals we used an antibonding 
combination of Gaussians with a r.m.s. width of 1 bohr. 
Each Gaussian was sitting halfway between one of the 
two atoms and the center of their common bond. Dur- 
ing the minimization fii decreased from 106.76 bohr^ to 
87.47 bohr^, having reached the minimum in less than 30 
steps. (An alternative is to choose the initial subspace 
at each k as the lowest four energy eigenstates inside the 
energy window. This yields an initial fli = 98.10 bohr^, 
and again the absolute minimum is reached after '^30 
steps.) The total spread of the four "maxloc" WFs was 
fl — 97.49 bohr^; as expected,L3 this is considerably 
larger than for the bonding WFs. Note also that fl ac- 
counts for more than 10% of the total spread, whereas for 
the bonding "maxloc" WFs that number was less than 
6%. This is related to the fact that the antibonding WFs 
are more spread out, causing matrix elements of the type 
(wmR|r|w„o) with R ^ to have larger values. Eq. (15) 
of Ref. |l6| shows that this results in a larger ft. The very 
small contribution of fl to the total spread of the highly 
localized d-like WFs in copper (less than 1%), as well as 
the comparatively larger contribution in the interstitial- 
centered WFs are thus easily understood. 
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FIG. 8. Solid lines: Original band structure of silicon. 
Dotted lines: Wannier-interpolated bands. In (a) the valence 
and low-lying conduction bands are treated separately, which 
produces four bonding and four antibonding Wannier func- 
tions; in (b) they are treated as a single group, which yields 
eight sp'^-type Wannier functions. 



In Fig. g(a) we present the contour-surface plot of one 
"maxloc" antibonding WF in silicon. The other three are 
identical (related to the first by the tetrahedral symmetry 
operations). Fig. ^(b) shows one of the four identical 
bonding WFs. 



2. sp hybrids 

As discussed in Ref. ^ one may instead treat the four 
valence and four low-lying conduction bands as a sin- 
gle group, which leads to "maxloc" WFs of sp^ char- 
acter (Fig. ^(c)). Using our method this may be done 
as indicated in Fig. ^(b). An outer energy window is 
chosen which spans the eight bands of interest, and the 
valence bands are "frozen" inside an inner window; this 
ensures that they are not affected by the minimization 
of Qi, whose only aim is to extract the four low- lying 
antibonding bands from the conduction band complex. 
We have started the minimization of Slj in two different 
ways: (i) by projecting eight "atom-centered" sp'^-type 
combinations of Gaussians, and (ii) by projecting four 
bond-centered Gaussians plus four antibonding combi- 
nations of Gaussians, as done in the previous Section. In 
both cases the minimization took about 20 steps, taking 
from 76.04 bohr-^ in the former case and 84.08 bohr'^ in 
the latter to 63.50 bohr^. As for the minimization of ft, 
the absolute minimum (il = 85.41 bohr^) was reached 
only with (i); with (ii) the algorithm became trapped in 
a local minimum (fi = 101.97 bohr^) having the same 
symmetry as the trial orbitals, with four bonding (anti- 
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FIG. 9. Contour-surface plots of Wannier functions in sil- 
icon, (a) Antibonding, (b) bonding, and (c) sp'^-type. In 
(a) and (c) the amplitudes are +Q.^/^/v (light gray) and 
—Q.^/^/v (dark gray); in (b) they are -1-1.4/y'u (light gray) 
and —QA/y/v (dark gray). 

bonding) WFs with a spread of 6.37 bohr^ (19.12 bohr^) 
each. 

We end this section with the following observation. 
Suppose we take the four-dimensional valence (bonding) 

space iS^'''* (k) together with the optimal four-dimensional 

antibonding subspace 54"' (k) (Fig. ||(a)) to form an 

eight-dimensional space <Sg (k) = 54^*^ (k) U (k) . This 
space has Vti = 63.64 bohr^, which is slightly higher than 
the value 63.50 bohr^ associated with the optimal sub- 
space iS8(k) for the eight-band problem with an inner 
window (Fig. ^(b)). Thus, if we take 58(k) as an initial 
guess for the minimization of fii in the eight-band prob- 
lem with an inner window, we will be starting slightly 
above the absolute minimum. The extra reduction in f2i 
comes about because the functional that is minimized to 
obtain tSslk) contains terms involving overlap between 
low-lying conduction states at k and valence states at 



neighboring k -I- b. The wavefunctions relax in response 
to these extra terms, and consequently the two antibond- 
ing subspaces are not exactly the same. However, they 
are almost identical, and therefore the same is true for 
the interpolated bands (compare Figs, ^(a) and||(b)). 

V. CONCLUSIONS 

We have discussed and implemented a practical 
method for extracting maximally-localized Wannier func- 
tions from entangled energy bands, starting from the 
Bloch eigenfunctions obtained in a standard electronic 
structure calculation. Our method is based on a pre- 
scription for "disentangling" the bands of interest from 
the rest of the band complex inside an energy window 
specified by the user. The idea is to extract a subspace 
of Bloch-like states whose character varies as little and 
as smoothly as possible across the Brillouin zone. This is 
achieved by minimizing a functional which measures the 
"spillage" , or change of character of the subspace across 
the Brillouin zone. The present scheme can be viewed as 
an extension of the maximally-localized-. Wannier func- 
tion method of Marzari and Vanderbilt,E3 which was de- 
signed to deal with isolated groups of bands only. More 
precisely, it introduces an extra step - the construction of 
the optimal subspace ~ which is followed by the determi- 
nation of the "maxloc" WFs by applying the localization 
algorithm of Marzari and Vanderbilt to that subspace. 
The procedure for determining this optimal subspace is 
both stable and computationally very fast. 

Some possible applications of such WFs have been 
mentioned in the Introduction. Of particular interest 
is the ability to obtain WFs for the low-lying empty or 
partially filled bands. For instance, it has been suggested 
that these could be useful for accurate calculations-af the 
optical properties of semiconducting nanocrystalsO An- 
other potential use of the present method could arise in 
the description of surface states (e.g., Ref. in par- 
ticular when the surface bands become resonant with 
the bulk bands. The striking result that we have ob- 
tained for the low- lying broad bands of copper, with the 
WFs being centered at the tetrahedral interstitial sites, 
suggests that the method may provide insight into the 
chemistry of transition metal compounds. Also, since the 
"maxloc" WFs provide a compact interpolation scheme 
for the band structure, they could be used as part of an 
efficient algorithm for determining the Fermi surface. Fi- 
nally, it might be interesting to apply the present ideas to 
the construction of lattice WFs describing the part of the 
phonon spectmm relevant for studying structural phase 
transitions .oEJ 
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